Distribution of Mechanical Properties in Poly(ethylene oxide)/silica Nanocomposites via Atomistic Simulations: From the Glassy to the Liquid State

Polymer nanocomposites exhibit a heterogeneous mechanical behavior that is strongly dependent on the interaction between the polymer matrix and the nanofiller. Here, we provide a detailed investigation of the mechanical response of model polymer nanocomposites under deformation, across a range of temperatures, from the glassy regime to the liquid one, via atomistic molecular dynamics simulations. We study the poly(ethylene oxide) matrix with silica nanoparticles (PEO/SiO2) as a model polymer nanocomposite system with attractive polymer/nanofiller interactions. Probing the properties of polymer chains at the molecular level reveals that the effective mass density of the matrix and interphase regions changes during deformation. This decrease in density is much more pronounced in the glassy state. We focus on factors that govern the mechanical response of PEO/SiO2 systems by investigating the distribution of the (local) mechanical properties, focusing on the polymer/nanofiller interphase and matrix regions. As expected when heating the system, a decrease in Young’s modulus is observed, accompanied by an increase in Poisson’s ratio. The observed differences regarding the rigidity between the interphase and the matrix region decrease as the temperature rises; at temperatures well above the glass-transition temperature, the rigidity of the interphase approaches the matrix one. To describe the nonlinear viscoelastic behavior of polymer chains, the elastic modulus of the PEO/SiO2 systems is further calculated as a function of the strain for the entire nanocomposite, as well as the interphase and matrix regions. The elastic modulus drops dramatically with increasing strain for both the matrix and the interphase, especially in the small-deformation regime. We also shed light on characteristic structural and dynamic attributes during deformation. Specifically, we examine the rearrangement behavior as well as the segmental and center-of-mass dynamics of polymer chains during deformation by probing the mobility of polymer chains in both axial and radial motions under deformation. The behavior of the polymer motion in the axial direction is dominated by the deformation, particularly at the interphase, whereas a more pronounced effect of the temperature is observed in the radial directions for both the interphase and matrix regions.


S1 Force field description for the PEO/SiO 2 model and atomistic simulation
For the description of the PEO, a modified united atom (UA) TraPPE based force field was used, 1,2 and for the SiO 2 nanoparticle a full atom representation was used. 3For the calculation of the electrostatic interactions the particle-mesh Ewald (PME) method was applied. 4The simulations were performed in the NPT statistical ensemble, where the pressure was kept constant with the use of Parrinello Rahman barostat 5 and the temperature is applied using the Nose Hoover thermostat. 6The LAMMPS simulation package was used for 1 the simulations. 7All force field parameters for PEO and silica are provided in in Tables S1   and S2.
Table S1: Non-Bonded Interactions (i, j denote different atom types)   To avoid "system size effects" in our simulations, i.e., the overall mechanical behavior and the properties computed (e.g.density profile, mechanical properties, chain dimensions, etc.) are the same if larger model systems, with the same loading of nanofillers, are used.To avoid such effects we considered large systems involving up to 8 nanoparticles.In order to improve statistics, all model atomistic PEO/SiO2 systems considered in the present work include 8 SiO2 nanoparticles and 384 PEO chains (each chain contains 150 atoms), whereas the size of the simulation cubic box at equilibrium (for one cell contains one NP) , in each direction, varies from 5.7nm at 220 K to 5.9 nm at 400 K.The counter length is 230 Å, Kuhn length is 10 Å, persistence length is 9.3 Å, 8 and radius of gyration in function of temperature are given in table S3.

S2 Thermodynamic equilibrium
To verify that the thermodynamic equilibrium is established, we monitored the potential energy, the density as well as the polymer radius of gyration.Figure S1 displays the corresponding graphs for a simulation run of 50 ns at a temperature T = 370K.Moreover, we perform additional short MD runs under NPT condition for few ns for thermal and local structure equilibration.To clarify this issue we present in Figure S2 the evolution of the energy (potential and kinetic energies) as well as the pressure during the relaxation for the replicated system at 220 K.The average values of the pressure are close to the atmospheric values.Large fluctuations in pressure are due to the incompressibility of the system and it is a common issue in MD simulations.

S3 Computation of local stress-strain
To investigate the heterogeneous mechanical behavior of the model PEO/SiO 2 systems we need to probe stress and strain field in local (atomic) level.For this, we use a per-atom calculation of stress and strain, under an imposed global strain.Local (per-atom) stress can be directly computed, for each atom i, σ i , via the atomic Virial formalism through: where the indices i, j, k and l denote atoms and α and β the Cartesian coordinate system (x, y, z).Eq. (S1) contains contributions from pairwise non-bonded interactions, bond stretching, bond angle bending, and dihedral torsional potentials, a K space contribution from long-range Columbic interactions and finally there is a term for the N f fixes that apply internal constraint forces to atom i. V i is the atomic volume of atom i.Based on the above expression, the stress within a specific domain (e.g.PEO/SiO 2 interface and PEO matrix) can be computed by summing up the per-atom contribution of all atoms within the domain and dividing with its volume.
The spatial distribution of strain fields at an atomic level can be probed directly for any system under study.In this approach, the strain distribution is computed by using concepts from continuum mechanics to define the strain per atom, relating the kinematics and the positions of atoms at an initial time t 0 and at a subsequent time t 1 .The coordinate vectors of atoms in the initial reference configuration Ω 0 at time t 0 having a Cartesian basis ẽi , and in the configuration Ω 0 at time t 1 with a Cartesian basis e i are represented by X and x , respectively ( Figure S3 ).Then the mapping from position vector X to x , is denoted by Figure S3: Kinematics of neighboring atoms at different times.A region in space is shown, where X and x represent the coordinate vectors of atoms in the reference configuration Ω 0 at time t 0 , and in the configuration Ω 1 at time t 1 , respectively.
for some coordinates X i and x i in their respective basis vectors.Assuming sufficient continuity, the local deformation at the spatial point X is characterized as the gradient of the mapping function motion and defined as where ∇ X denotes the gradient operator taken in with respect to X.The deformation gradient F obtained from Eq. ( S3) is applicable at atomic scales and is also available for the continuum framework.
It is expected that the deformation in the neighborhood of atom m is characterized by the changes in the relative position of its neighbor n.Atom m is located at position X m in the reference configuration Ω 0 and position x m in the current configuration Ω 1 .Then, the relative position of neighboring atom n and the deformation gradient F m at atom m is related to its neighboring atoms, and satisfies where • denotes dot product tensor contractions and In principle, for a given reference atom m, Eq. (S4a) must hold for all its neighboring atoms.
However, F m of atom m cannot be generally determined from a single atom n, since, unlike the continuum limit, the interparticle interactions with all its neighboring atoms need to be taken into account, as their motion is strongly correlated.Therefore, we seek an optimal local deformation gradient Fm (x), which is obtained from the minimization problem: i.e. it minimizes the squared error where N is the number of neighboring atoms within the cutoff radius r cut and ∥ • ∥ denotes the Euclidean norm.From the optimal deformation gradient matrix Fm (x), we can then define the Lagrange Green strain tensor ε m (x) with respect to the reference coordinates as: where, for each atom m within the simulation box, ε m is a symmetric 3 × 3 tensor.

S4 Structure and Dynamics of Polymer Chains in the Interfacial Zone
Regarding the influence of silica nanoparticles on polymer dynamics, in a recent work by our group, we examined the spatial heterogeneous dynamics and glass transition of polymer chains in PEO/SiO 2 nanocomposites with varying silica concentration. 9The results demonstrate the existence of a spatial gradient in the segmental dynamics and the glass transition temperature, as well.In more detail, a layer-resolved analysis of the mean squared displacement (MSD) of the polymer chains revealed that the polymer chains belonging to the matrix region are more mobile than their counterparts residing in the interphase region as it can be seen in Figure S4 where the MSD of the interphase and matrix regions is plotted as a function of time at temperatures T = 300 and T = 400K.
Focusing on Figure S4, we can also deduce that the slowing down in polyner mobility becomes more pronounced upon approaching the glassy regime (T = 300K).This implies that confinement effects have a higher impact on the polymer dynamics in the glassy state.

S5 Determination of glass transition temperature
The transition temperature is a critical property in our study.Actually, the glass transition temperature of the examined system has already been calculated recently by our group. 9

Figure S1 :Figure S2 :
Figure S1: Time evolution of the potential energy (left), density (right) and radius of gyration (bottom) of the system T= 370 K.

Figure S4 :
Figure S4: Mean squared displacement (MSD) of the polymer chains residing at the interphase or matrix regions at T= 300, 400 K. .

FigureFigure S5 : 7 S p e c i f i c V o l u m e S p e c i f i c v o l u m e ( c m 3 /Figure S6 :
Figure S6 illustrates the temperature dependence of the specific volume of the nanocomposite

Table S3 :
Evolution of R g in function of temperature